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Abstract 
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the  time-dependent,  three-dimensional,  arbitrary-geometry  MHD  model  which  includes 
viscous  and  resistive  effects.  A  time-dependent  ionization  model  was  added  which  self- 
consistently  calculates  the  ionization  fraction  of  the  fluid.  Energy  loss  mechanisms 
were  added  for  the  constituent  fluid  components  (neutrals,  ions,  and  electrons).  The 
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1  Executive  Summary 

The  primary  objective  of  this  project  is  to  develop  a  novel  algorithm  to  accurately  simu¬ 
late  realistic  plasmas  with  constituent  species  at  multiple  temperatures  using  the  MHD 
model.  A  viable  time-dependent,  three-dimensional  MHD  code  will  provide  a  valuable 
tool  for  the  design  and  testing  of  plasma  related  technologies  that  are  important  to  the 
Air  Force  and  industry,  such  as  portable  pulsed  power,  high  power  microwave  devices, 
hypersonic  drag  reduction,  advanced  plasma  thrusters  for  space  propulsion,  nuclear 
weapons  effects  simulations,  radiation  production  for  counter  proliferation,  and  fusion 
for  power  generation.  Implementing  the  algorithm  on  parallel  supercomputers  will  allow 
the  detailed  modeling  of  realistic  plasmas  in  complex  three-dimensional  geometries. 

Current  MHD  codes  are  limited  to  simulations  of  short  time  scale  phenomena  be¬ 
cause  of  explicit  time  step  stability  limitations  and  equation  decoupling.  We  devel- 
opeedd  an  implicit  algorithm  with  the  capability  to  simulate  physics  of  any  length  time 
scale  because  the  time  step  is  chosen  by  the  user  to  match  the  physics  of  interest.  This 
algorithm  has  the  additional  advantage  that  the  equations  are  solved  in  a  fully  coupled 
manner.  The  plasma  is  assumed  to  be  composed  of  a  neutral  fluid,  ion  fluid,  and  elec¬ 
tron  fluid.  Each  fluid  has  an  associated  temperature  and  can  exchange  energy  to  the 
other  fluids  by  ionization  and  other  collisional  processes.  The  plasma  is  allowed  to  have 
a  variable  degree  of  ionization,  from  a  fully  ionized  plasma  to  a  completely  neutral  gas. 
The  algorithm  is  implemented  using  arbitrary  finite  volumes  so  it  can  model  realistic 
three-dimensional  geometries. 

To  speed  development  of  this  effort  an  existing  MHD  code  was  used.  WARPS  (Wash¬ 
ington  Approximate  Riemann  Plasma  code  for  3-d  domains)  is  a  time-dependent,  three- 
dimensional,  arbitrary-geometry  MHD  algorithm  with  viscous  and  resistive  effects.  The 
code  was  extended  to  include  thermal  diffusion  for  the  constituent  temperatures  (neu¬ 
trals,  ions,  and  electrons).  A  time-dependent  ionization  model  was  added  which  self- 
consistently  calculates  the  ionization  fraction  of  the  fluid.  Energy  loss  mechanisms  were 
added  for  the  constituent  fluid  components.  These  features  were  benchmarked  against 
analytical  results.  The  new  algorithm  is  solved  using  a  domain  decomposition  technique 
for  parallel  computation. 

The  implicit  formulation  has  been  developed  for  the  resistive  and  viscous  MHD 
model.  The  culmination  of  this  research  effort  produced  the  Ph.  D.  dissertation  of  B. 
Udreafl]  and  the  M.  S.  thesis  for  W.  Vuillemot.  The  algorithm  has  been  cast  using 
finite  volumes  which  significantly  reduces  transverse  flux  errors.  An  important  result 
of  this  work  is  the  development  of  a  2-level  nested  iteration  technique  which  accurately 
solves  the  MHD  equations  with  typical  Courant  numbers  of  100.  The  residual  of  the 
error  is  driven  to  machine  accuracy  for  all  cases  investigated. 

As  a  result  of  this  project  several  professional  collaborations  now  exist  between  the 
Department  of  Aeronautics  and  Astronautics  at  the  University  of  Washington  and  the 
Air  Force  Research  Laboratory,  Lawrence  Livermore  National  Laboratory,  the  Univer¬ 
sity  of  Michigan,  the  University  of  Colorado,  Stanford  University,  and  other  depart¬ 
ments  at  the  University  of  Washington.  The  work  from  this  project  has  been  presented 
at  international  conferences  and  published  in  a  refereed  journal. 


2  Project  Description 

Plasmas  are  essential  to  many  technologies  that  are  important  to  the  Air  Force,  some 
of  which  have  dual-use  potential.  These  applications  include  portable  pulsed  power  sys- 
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terns,  high  power  microwave  devices,  drag  reduction  for  hypersonic  vehicles,  advanced 
plasma  thrusters  for  space  propulsion,  nuclear  weapons  effects  simulations,  radiation 
production  for  counter  proliferation,  and  fusion  for  power  generation.  Several  of  these 
applications  are  specifically  mentioned  in  the  New  World  Vistas  Report  from  the  USAF 
Scientific  Advisory  Board. [2]  In  general,  plasmas  fall  into  a  density  regime  where  they 
exhibit  both  collective  (fluid)  behavior  and  individual  (particle)  behavior.  Many  plas¬ 
mas  of  interest  can  be  modeled  by  treating  the  plasma  like  a  conducting  fluid  and 
assigning  macroscopic  parameters  that  accurately  describe  its  particle-like  interactions. 
Magnetohydrodynamic  models  the  plasma  in  this  manner. 


2.1  Research  Objectives 

The  objectives  of  the  project  are  to: 

•  Develop  an  implicit,  conservative  multi-temperature  algorithm  for  three-dimensional 
non-ideal  MHD  simulations  for  time-dependent  and  steady  state  variably  ionized 
plasmas; 

•  Validate  the  code  with  analytical  and  experimental  data;  and 

•  Apply  the  code  to  analyze  plasma  related  topics  at  the  Air  Force  Research  Lab¬ 
oratories  [the  magnetic  flux  compression  generator  (MCG)  experiments  and  the 
liner  implosion  system  (WFX)[3]  at  Kirtland  AFB,  the  plasma  thruster  work  at 
Edwards  AFB,  and  the  hypersonic  drag  reduction  research  at  Wright-Patterson 
AFB]  and  at  the  University  of  Washington  [Z-Pinch  experiment  (ZaP)[4]  and  He- 
licity  Injected  Tokamak  (HIT)  [5]]. 


2.2  Technical  Description 

The  three-dimensional,  extended  MHD  plasma  model  is  a  set  of  mixed  hyperbolic  and 
parabolic  equations.  The  Navier-Stokes  equations  are  also  of  this  type.  This  project 
applies  some  advances  that  have  been  made  in  implicit  algorithms  for  the  Navier-Stokes 
equations  to  the  MHD  equations.  These  implicit  algorithms  solve  the  equation  set  in  a 
fully  coupled  manner,  which  generates  better  accuracy  than  the  current  methods  used 
for  MHD  simulations. 

When  expressed  in  conservative,  non-dimensional  form,  the  MHD  model  is  described 
by  the  following  equation  set. 
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The  variables  are  density  (p),  velocity  (v),  magnetic  induction  (B),  ionization  fraction 
(/i)>  pressure  (p),  and  energy  density  (e).  The  energy  density  is 
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where  7  =  cp/cy  is  the  ratio  of  the  specific  heats.  The  pressure  is  the  total  material 
pressure,  which  is  the  sum  of  the  partial  pressures  from  the  neutral,  ion,  and  electron 
fluids. 

p  =  nnkTn  +  riikTi  +  nekTe  (3) 

where  k  is  the  Boltzmann  constant,  nn  is  the  neutral  number  density,  and  Tn  is  the 
neutral  temperature.  The  remaining  variables  are  the  ion  and  electron  number  densities 
and  temperatures. 


p  =  (nn  +  Ui)Mi 


(4) 


and 


Tie  -  Tli. 


The  number  densities  are  determined  from  a  time-dependent  ionization  model 
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where  <  av  >ion  is  the  ionization  rate  parameter  and  <  av  > recomb  is  the  recombina¬ 
tion  rate  parameter.  [6]  The  multiple  temperatures  evolve  independently  based  on  the 
appropriate  components  of  the  energy  equation  and  energy  transfer  between  species. 

The  right  hand  side  of  eqn(l)  contains  the  non-ideal  effects.  These  effects  include 
viscosity,  resistivity,  Hall  currents,  diamagnetic  currents,  thermal  conduction,  and  ra r 
diation  cooling.  The  non-ideal  terms  are  defined  by 


PVvisc  = 


ReAl 


V-f 


(7) 


Bre.  =  - 


1 


RmAl 


Vx(?.(Vx  B)) 


—  WceTe  T7  v 

Bhm  =  -  w/v  x 


(V  X  B)  X  B 


(8) 

(9) 


Pvisc  =  ,,V  *  V  *  T 


Prea  = 


1 


RmAl 


ReAl 


V  •  tj  ■  (V  x  B)  x  B 


(10) 

(U) 

(12) 


Pi cond 


Mi 

2PeAl 


V  •  k  ■  (VTn  +  VTt  +  VTe)  =  Pcondn  +  P condi  +  Pcond.  (13) 


3 


Prad  =  -CradZ'B(pft)2TlJ2  (14) 

Mi  is  the  ion  mass,  ojcere  is  the  Hall  parameter,  Crad  Is  the  Bremsstrahlung  radiation 
constant,  and  Zeff  is  the  effective  ionization  level  due  to  plasma  impurities. 

The  non-dimensional  tensors  are  the_  stress  tensor  (f),  the  electrical  resistivity  (^), 
and  the  thermal  conductivity  (fc),  and  I  is  the  identity  matrix.  The  non-dimensional 
numbers  are  defined  as  follows: 

Alfven  Number  : 

Reynolds  Number  : 

Magnetic  Reynolds  Number 
Peclet  Number  : 

The  characteristic  variables  are  length  (L),  velocity  (V),  Alfven  speed  (Va  =  B/ y/p0p), 
kinematic  viscosity  (*✓),  electrical  resistivity  (rj),  and  thermal  diffusivity  (k  =  k/pCp ). 
fi0  is  the  permeability  of  free  space  (47r  x  10“7). 

For  convenience,  the  MHD  equation  set  [eqn(l)]  is  rewritten  in  the  following  compact 
form 

^  +  V  ■  f  h  =  Q  Non-ideal  (16) 

where  Q  is  the  vector  of  conservative  variables,  is  the  tensor  of  hyperbolic  fluxes, 
and  Q Non- ideal  contains  the  non-ideal  (mostly  parabolic)  terms.  The  forms  of  these 
vectors  and  tensors  can  be  seen  from  the  previous  equations.  The  hyperbolic  fluxes  are 
associated  with  wave-like  motion,  and  the  parabolic  fluxes  axe  associated  with  diffusion¬ 
like  motion. 


2.3  Multiple  Temperature  Evolution 

The  temperatures  of  the  multiple  species  (neutrals,  ions,  electrons)  evolve  indepen¬ 
dently  based  on  the  appropriate  components  of  the  energy  equation  and  energy  transfer 
between  species.  The  temperatures  are  consistent  with  energy  conservation. 

The  temperature  rise  in  each  specie  depends  on  the  heating  mechanism  and  the 
density  fraction  of  the  specie.  We  define  the  density  fractions  as 


fi  = 


m 

4*  Tii 


(17) 


(18) 

(19) 


For  multiply  charged  species,  the  last  definition  would  be  modified  to  fe  =  Zfc.  Viscous 
drag  heats  the  neutral  gas  and  the  ion  fluid,  but  does  not  affect  the  electrons.  Therefore, 
the  energy  rise  due  to  viscous  effects  is  attributed  to  the  neutral  and  ion  temperatures. 
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Resistivity  directly  heats  only  the  electrons. 

=  (7  —  l)^rej 
Radiation  is  emitted  by  the  electrons  as  they  cool. 

=  (7  -  l)Prad 
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Each  species  has  its  own  thermal  conduction  component  as  is  evident  from  eqn(13). 
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The  remaining  total  material  pressure  rise  is  due  to  adiabatic  compression  which  affects 
all  species  in  proportion  to  their  densities. 

Energy  may  also  transfer  between  species  due  to  collisions.  However,  total  energy 
is  conserved  and  the  total  material  pressure  is  not  affected.  The  interspecies  energy 
transfer  is  modeled  to  preserve  the  total  material  pressure. 


P  =  Pn+Pi  +  Pe  =  Pn+Pi  +  Pe 


(27) 


2.4  Conservative  Algorithm 

Because  of  the  natural  differences  between  hyperbolic  and  parabolic  equations,  the 
methods  for  solving  them  are  very  different.  Since  the  MHD  equations  are  of  mixed 
type  the  hyperbolic  and  parabolic  terms  must  be  handled  differently.  The  hyperbolic 
fluxes  are  differenced  by  applying  an  implicit,  approximate  Riemann  algorithm  that 
properly  accounts  for  their  wave-like  behavior.  The  parabolic  terms  are  discretized  by 
applying  explicit  central  differencing.  The  remaining  non-ideal  terms  which  correspond 
to  the  Hall  effect  are  solved  using  a  semi-implicit  method.  [7] 

The  design  of  the  algorithm  is  driven  by  the  conservative  numerical  techniques  that 
must  be  used  for  the  hyperbolic  terms.  Therefore,  we  begin  by  considering  the  ideal 
MHD  equations,  which  are  obtained  from  eqn(16)  by  setting  all  the  non-ideal  terms 
(Q Non— ideal)  to  zero.  Note  that  ideal  MHD  refers  to  an  ideal  plasma  —  one  that  is 
inviscid  and  non-resistive  and  neglects  thermal  conduction  and  finite  Larmor  radius 
(FLR)  effects. 

The  ideal  MHD  equations  are 

^  +  V  '  +  /»V  ■  <5  =  0,  (28) 
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where  A  is  the  Jacobian  of  the  hyperbolic  flux  tensor. 


Here,  Q  is  the  vector  of  conserved  variables. 

Q  =  (pypvx,pVy,pvz,BXiByyBZ)e)  . 


(29) 


(30) 


This  is  a  set  of  hyperbolic  equations  and  thus  Ax  has  a  complete  set  of  real  eigenvalues 
given  by 

A  —  (Vx>0,l>£  ±  V/asti^x  i  VslowjVx  i  Vi4x)  j  (31) 


where  Vfast  and  V3iow  are  the  fast  and  slow  magnetosonic  speeds,  and  Vax  is  the  Alfven 
speed  based  on  the  x  component  of  the  magnetic  field.  These  can  be  expressed  as 
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c*  +  V t+  \J {c2s  +  Vlf  -  4c2sV%x 


K2 


slow 


c2  +  K,2  -  yf{&- 7vff-^v[x 


(32) 

(33) 


Y  Ax 


Bl 

PoP' 


Here,  cs  is  the  ion  sound  speed,  which  for  a  perfect  gas  is 


7 P 
P  ' 


(34) 


(35) 


We  make  special  note  of  the  zero  eigenvalue,  A2  in  this  case.  The  zero  eigenvalue 
only  appears  in  multiple  dimensions  and  is  caused  by  the  perpendicular  nature  of  the 
j  x  B  force.  Powell  et  al\ 8]  recently  solved  this  zero  eigenvalue  problem  by  introducing 
a  source  term  that  is  proportional  the  divergence  of  the  magnetic  field.  The  eigenvalue 
becomes  finite,  A2  —  vx  in  this  case.  We  have  implemented  this  modification  and  it  is 
discussed  in  a  later  section. 

For  hyperbolic  equations  information  propagates  along  characteristics  which  travel 
at  wave  speeds  given  by  the  eigenvalues.  Most  modern  numerical  techniques  for  solving 
hyperbolic  equations  are  based  upon  the  idea  of  splitting  the  fluxes  into  components 
due  to  left-  and  right-running  waves.  Then  each  part  of  the  flux  can  be  differenced 
in  an  upwind  manner,  which  greatly  reduces  numerical  oscillations  and  stabilizes  the 
solutions. 

It  is  well  known  that  if  a  hyperbolic  equation  is  solved  with  an  explicit  scheme,  then 
the  allowable  time  step  to  maintain  numerical  stability  is  given  by  the  CFL  (Courant- 
Friedrichs-Lewy)  condition,  which  in  the  case  of  the  ID  MHD  equations  is 


A  t  < 


Ax 

I^X  j-  Vfast\ 


(36) 


For  the  high  magnetic  fields  and  low  densities  common  in  many  plasma  experiments,  the 
fast  magnetosonic  speed  is  quite  high,  and  thus  the  time  step  is  prohibitively  small.  We 
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are  often  interested  in  only  modeling  the  physics  that  occurs  slower  than  Alfven  time 
scales.  For  example,  it  can  be  shown  that  resistive  tearing  modes,  which  are  important 
in  studying  fusion  plasmas,  evolve  on  a  time  scale  that  is  given  by [9] 

Ttearing  OC  =  (Luf/5  TA.  (37) 

ta  is  the  Alfven  time,  rn  is  the  resistive  diffusion  time,  and  Lu  is  the  Lundquist  number, 
which  is  given  by 

Lu  =  &  =  RmAL  (38) 

ta 

If  Lu  is  106,  which  is  typical  for  laboratory  plasmas  in  fusion  applications,  the  resistive 
tearing  time  is  approximately  4000  times  larger  than  the  Alfven  time.  By  treating  the 
hyperbolic  fluxes  implicitly  in  time,  the  stability  restriction  on  the  time  step  is  removed, 
and  the  solution  can  be  advanced  at  the  larger  resistive  tearing  time  step.  This  is  our 
motivation  for  proposing  an  implicit  scheme. 


2.5  Implicit  Formulation 

For  clarity,  the  algorithm  for  the  two-dimensional  ideal  MHD  equations  is  presented. 
The  extension  to  three  dimensions  is  straight  forward.  The  two-dimensional  equation 
is 


?Q  +  dF  +  dG=zQ 

dt  dx  dy 


(39) 


Eqn(39)  was  discretized  using  first  order  Euler  time  differencing  to  get 

(Q£+ 1  -  Q?jf)  -  -Rn  (Qn+1)  =  - R\ 


(40) 


where  R  is 


Rij  =  Fi+i/zj  ~  Ft- 1/2, j  +  Gij+1/2  -  Gij- 1/2.  (41) 

Note  that  in  this  equation  and  all  that  follow  the  grid  metric  terms  (cell  areas  and 
volumes)  were  omitted  for  clarity.  Linearizing  R  as  follows: 

K+l «  Rij  +  (^)".  (%+1  -  QW  (42> 


where  dR/dQ  has  been  defined  as  the  differenced  flux  Jacobians. 


dR 

dQ 


dF 

dQ 


i+l/2  ,j 


dF 

dQ 


+  ... 


(43) 


where  dF/dQ  is  the  flux  Jacobians  of  the  x  flux.  The  flux  Jacobians  can  be  calculated 
analytically  or  numerically.  Analytical  calculations  based  on  the  assumption  that  the 
solution  values  do  not  change  rapidly  were  used  previously  [10]  and  produced  adequate 
results.  The  current  project  investigated  numerical  calculation  of  the  flux  Jacobians  and 
an  analytical  method  without  the  previous  assumption.  Two  methods  for  numerical 
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calculation  were  investigated.  First,  a  limit  formulation  similar  to  the  definition  of  a 
differential  was  used. 


dF  F{Q  +  e)-F(Q)  () 
dQ  e  K> 


(44) 


for  small  e  which  gave  first  order  accuracy  in  e.  The  flux  Jacobians  were  also  calculated 
using  complex  numbers.  The  flux  Jacobians  were  expanded  about  Q  in  a  Taylor  series. 


The  expression  was  rearranged  to  solve  for  the  flux  Jacobian. 


dF 

dQ 


F(Q  +  ih) 
h 


+  0{h2) 


(46) 


which  gave  second  order  accuracy  in  h.  Additionally,  the  complex  formulation  required 
only  a  single  evaluation  of  the  flux  Jacobian  (though  using  complex  math)  compared  to 
two  evaluations  for  the  limit  formulation. 

Substituting  the  expression  for  R”+l  back  into  eqn(40)  and  rearranging,  gave 


A  Q"  =  ~R?y 


(47) 


Here  AQ  has  been  defined  as 

A  Qn  =  Q"/1  -  QS- 


(48) 


The  left  hand  side  of  the  eqn(47)  is  an  implicit  operator  operating  on  A Q.  It  is  a 
large  banded  block  matrix.  In  three  dimensions,  it  is  an  {I max  X  Jmax  x  Kmax)  hy 
(Imax  x  Jmax  x  Kma x)  matrix,  where  Imax  is  the  number  of  cells  in  the  x  direction,  etc. 
It  is  quite  costly  to  invert  a  matrix  of  this  size  directly.  For  this  project  a  more  efficient 
iterative  method  was  chosen.  When  solved  iteratively  ,  eqn(47)  can  lose  time  accuracy. 
To  recover  time  accuracy  the  time  derivative  of  Q  was  added  as  a  source  term  to  the 
right  hand  side  of  the  equation.  The  modified  equation  became 


—  — .RTV*"1  _  gn+l 


(49) 


where 

S"+1  =  2^  (3Q"+1  -  4Q<j  +  Qif1)  *  f^.  (50) 

The  r  in  eqn(49)  is  an  iteration  variable, called  pseudo  time.  At  each  physical  time 
step,  eqn(49)  is  solved  iteratively  until  the  left  hand  side  vanishes.  When  the  solution 
converges,  our  original  equation 


dQ 

dt 


=  -R 


(51) 


is  solved.  This  technique  is  known  as  dual  time-stepping.[ll]  Note  that  in  eqn(50)  a  more 
accurate  time  derivative  can  be  used  at  the  expense  of  the  additional  memory  needed 
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to  store  the  Q  vectors  from  previous  time  steps.  A  concern  for  higher  order  numerical 
methods  is  the  property  of  total  variation  diminishing  (TVD).  The  TVD  property  pre¬ 
vents  some  types  of  nonlinear  numerical  instabilities.  The  dual  time-stepping  method 
has  not  been  shown  to  be  TVD  for  the  MHD  model. 

One  advantage  of  the  strategy  outlined  above  is  that  the  implicit  operator  and  the 
right  hand  side  in  eqn(47)  are  decoupled.  The  structure  of  the  matrix  no  longer  depends 
on  the  details  of  the  discretization  of  the  right  hand  side  fluxes.  The  iteration  equation 
[eqn(47)] 

Sa^2a;  +  (ID"] =  -5S  + -  *3  (52) 

has  the  standard  form 


Ax  =  b.  (53) 

Previously  the  LU-SGS  method  (lower-upper  symmetric  Gauss-Seidel)[12]  was  used 
to  invert  the  implicit  operator  A.  [13,  10]  The  LU-SGS  method  required  a  modification 
of  the  implicit  operator  through  an  approximate  factorization  procedure  which  reduced 
the  accuracy  of  the  operator  and  led  to  poor  convergence. 

The  current  project  used  a  symmetric  Gauss-Seidel  method  which  does  not  rely  on 
approximate  factorization.  The  SGS  method  was  used  to  iteratively  invert  the  implicit 
operator  and  the  approximate  Riemann  solver  that  is  used  to  form  the  right  hand  side 
fluxes. 

The  implicit  operator  matrix  was  decomposed  into  lower,  upper,  and  diagonal  ma¬ 
trices  and  written  as 


A  =  L  +  U  +  D  (54) 

Each  iteration  of  the  symmetric  Gauss-Seidel  method  performs  two  sweeps  of  the  domain 
—  a  forward  sweep  followed  by  a  backward  sweep. 


(L  +  DJx^  +  Ux^-^b 

(55) 

(U  +  D)x(2I)  +  Lx(2'-1}  =  b 

(56) 

where  l  —  1,2,3,...  is  the  iteration  index. 

For  reasonably  large  values  of  Re  and  Rra  (easily  in  the  range  of  interest  for  most 
applications),  the  parabolic  terms  can  be  differenced  explicitly  without  constraining  the 
allowable  time  step.  The  right  hand  side  of  eqn(47)  was  modified  by  adding  the  central 
differenced  parabolic  terms. 

2.6  Approximate  Riemann  Solver 

The  fluxes  on  the  right  hand  side  of  eqn(47)  were  discretized  using  a  Roe-type  ap¬ 
proximate  Riemann  solver.  [14]  In  this  method  the  overall  solution  was  built  upon  the 
solutions  to  the  Riemann  problem  defined  by  the  discontinuous  jump  in  the  solution 
between  each  pair  of  cells.  The  numerical  flux  for  a  first-order  accurate  (in  space) 
Roe-type  solver  was  written  in  symmetric  form  as 

Fi+ 1/2  =  \  (Fi+x  +  F,)  -  I  £>  (Qi+I  -  Qi)  lAfel  rk  (57) 
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where  rk  is  the  kth  right  eigenvector,  A*  is  the  absolute  value  of  the  kth  eigenvalue,  and 
Ik  is  the  kth  left  eigenvector.  The  values  at  the  cell  interface  (i  4-  1/2)  were  obtained 
by  either  a  simple  average  or,  more  accurately,  a  Roe-average  of  the  neighboring  cells. 
Determining  a  Roe-average  on  an  arbitrary  computational  grid  involved  transforming 
the  vector  quantities  to  a  coordinate  system  that  is  orthogonal  to  the  local  cell  interface. 
Then  the  flux  calculated  as  above  is  normal  to  the  cell  interface  which  is  the  desired 
orientation  for  applying  the  divergence  theorem  in  a  finite  volume  method. 

These  first  order  accurate  upwind  fluxes  are  used  in  the  vicinity  of  sharp  disconti¬ 
nuities  in  order  to  suppress  oscillations  in  the  solution.  Globally  second  order  accurate 
solutions  were  achieved  by  using  a  flux  limiter  that  modifies  the  first  order  flux  so  that 
it  uses  second  order  central  differencing  in  smooth  portions  of  the  flow.  The  minmod 
limiter  was  used.  [15] 

Once  the  eigenvalues  and  eigenvectors  were  obtained  and  properly  normalized  to 
avoid  singularities,  it  was  relatively  straight-forward  to  apply  this  scheme  to  the  one¬ 
dimensional  ideal  MHD  equations.  [16,  17]  Unlike  for  the  Euler  equations,  the  extension 
to  more  than  one  dimension  was  non- trivial.  In  more  than  one  dimension,  the  Q  vector 
must  include  Bx  in  addition  to  the  other  magnetic  field  components.  (For  the  one¬ 
dimensional  case  Bx  is  constant  by  virtue  of  V  •  B  =  0).  Since  the  j  x  B  force  acts 
perpendicularly  to  the  directions  of  j  and  B,  the  F  flux  vector  has  a  zero  term  corre¬ 
sponding  to  Bx .  Thus,  the  Jacobian  matrix  of  F  is  singular  and  has  a  zero  eigenvalue. 
A  complete  set  of  physically  meaningful  eigenvectors  no  longer  exists.  Physically,  one 
would  expect  information  to  travel  either  at  the  fluid  velocity  or  at  the  fluid  velocity 
plus  or  minus  the  wave  speeds.  Simply  dropping  Bx  from  the  equation  set  is  not  a 
viable  option,  because  Bx  needs  to  change  in  order  to  maintain  V  •  B  =  0.  Powell  et 
af.,[8]  solved  this  problem  by  modifying  the  Jacobian  in  such  a  way  as  to  change  the 
zero  eigenvalue  to  vx  (keeping  the  others  unchanged) ,  and  then  adding  in  a  source  term 
that  exactly  canceled  out  the  terms  introduced  by  the  modified  Jacobian. 

The  source  term  is 


P 

B 

v 

V  B 


V  B 


(58) 


It  is  proportional  to  the  divergence  of  B  and  thus  very  small.  The  effect  of  this  mod¬ 
ification  is  to  sweep  the  field  divergence  out  of  the  domain  with  the  plasma  flow  (for 
example,  Afc  =  vx).  However,  for  closed  boundaries  or  stagnation  points  the  divergence 
increases.  A  divergence  cleaner  has  been  implemented  based  on  the  Hodge  projection 
technique. 

2.7  Finite  Volume  Grid  and  Parallel  Implementation 

Since  the  algorithm  being  developed  will  be  used  for  real  devices,  it  must  be  capable  of 
modeling  arbitrary,  three-dimensional  geometries.  Therefore,  multi-block,  finite  volume 
grids  were  used.  A  typical  cell  is  shown  in  Figure  1.  The  computational  domain  is 
divided  into  blocks  which  are  then  spanned  by  body-fitting  finite  volume  cells.  See 
Figure  2  for  a  possible  grid. 

As  discussed  in  the  previous  section,  the  formulation  of  the  approximate  Riemann 
solver  which  we  have  developed  generates  hyperbolic  fluxes  oriented  normal  to  the  grid 
cell  interfaces.  Application  of  the  divergence  theorem  is  then  a  simple  operation.  The 
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Figure  1:  Schematic  of  the  arbitrary  shaped  three-dimensional  finite  volume  cell  used  by 
the  algorithm. 


Figure  2:  Schematic  of  the  arbitrary  shaped  three-dimensional  finite  volume  cell  used  by 
the  algorithm. 
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parabolic  fluxes  are  also  calculated  to  be  normal  to  the  grid  cell  interfaces.  To  accom¬ 
plish  this  orientation,  a  set  of  nested  control  volumes  were  used  and  the  appropriate 
vector  operations  within  these  volumes  were  applied. 

The  block  structure  of  the  grid  provided  a  natural  domain  decomposition  for  the 
parallel  implementation.  The  integral  form  of  a  general  conservation  law  was  expressed 
as 

2-fdVQ+idS-  F(Q)  =  JdV  S(Q),  (59) 

ns  n 

where  ft  is  the  domain  and  E  is  the  boundary  of  ft.  Q  is  the  vector  of  conserved 
variables,  F(Q)  is  the  flux  of  the  conserved  variables,  and  S(Q)  is  the  vector  of  source 
terms.  By  splitting  the  domain  ft  into  p  subdomains  such  that 

n  =  U  Qu  (60) 

i  —  1 

eqn(59)  was  replaced  with  a  set  of  p  conservation  equations  applied  on  the  subdomains 

fti. 

jtJ  dV  Q  +  jdS-F(Q)  =  J  dV  S(Q),  i  =  l,2,...,p  (61) 

Ei  ft* 

Each  of  these  discretized  equations  is  solved  by  a  single  processor  using  the  boundary 
values  copied  from  neighboring  subdomains. 

To  ensure  a  portable  code  a  message  passing  system  commonly  available  on  parallel 
supercomputers  and  on  workstation  clusters  was  used.  This  system  was  the  Message 
Passing  Interface  (MPI)[18],  which  was  adopted  as  a  standard  in  May  1994  by  indus¬ 
try  and  academia.  Hardware  and  software  vendors’  implementation  of  MPI  provides 
parallel  program  developers  with  a  consistent  set  of  subroutines  callable  from  FOR- 
TRAN90  and  C.  In  this  project  the  basic  point-to-point  communications  subroutines 
and  global  communications  subroutines  were  used.  The  point-to-point  communication 
subroutines  were  used  for  the  domain  decomposition  and  boundary  exchange  while  the 
global  communication  subroutines  were  used  for  convergence  checking.  All  message 
passing  systems  (PVM,  MPL)  support  point-to-point  and  global  communications  sub¬ 
routines  so  that  by  using  only  the  basic  set  portability  to  systems  not  supporting  MPI 
was  simplified. 

3  Project  Implementation  and  Results 

3.1  Finite  Volume  Improvements 

To  improve  the  codes  ability  to  handle  highly  distorted  grids,  finite  volume  grids  were 
implemented.  The  finite  volume  implementation  greatly  reduced  the  anomalous  mo¬ 
mentum  leakage  into  orthogonal  directions  when  the  grid  was  distorted. 

Figure  3  shows  a  shock  tube  test  problem.  The  simulation  was  performed  in  three- 
dimensions,  but  should  remain  one-dimensional.  The  figure  shows  the  three  gas  dynamic 
waves  in  the  density  plot.  The  transverse  velocity  should  be  zero.  A  finite  amount  of 
momentum  leakage  was  generated  by  the  grid  metrics  in  the  finite  difference  generalized 
coordinate  formulation. 
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Figure  3:  A  three-dimensional  simulation  of  a  one-dimensional  shock  tube  showing  the 
density  and  transverse  velocity.  The  three  gas  dynamic  waves  can  be  seen  in  the  density 
plot.  The  transverse  velocity  should  be  zero. 


13 


Figure  4:  A  three-dimensional  simulation  of  a  one-dimensional  shock  tube  showing  the 
density  and  transverse  velocity.  The  three  gas  dynamic  waves  can  be  seen  in  the  density 
plot.  The  transverse  velocity  should  be  zero. 


The  reduction  of  the  anomalous  momentum  leakage  into  orthogonal  directions  can  be 
seen  in  Figure  4.  The  simulation  was  identical  to  that  shown  in  Figure  3  except  a  finite 
volume  implementation  was  used  instead  of  the  generalized  coordinate  formulation. 
Figure  4  shows  the  same  three  gas  dynamic  waves  in  the  density  plot.  With  the  finite 
volume  implementation,  the  transverse  velocity  is  zero  to  within  machine  accuracy. 

3.2  Implicit  Formulation  and  Numerical  Flux  Jacobian  Calcu¬ 
lations 

The  LU-SGS  method  (lower-upper  symmetric  Gauss-Seidel)  [12]  was  used  previously  to 
invert  the  implicit  operator  of  eqn(47).[13,  10]  The  LU-SGS  method  required  a  modifi¬ 
cation  of  the  implicit  operator  through  an  approximate  factorization  procedure  which 
reduced  the  accuracy  of  the  operator  and  led  to  poor  convergence.  The  convergence 
history  is  shown  in  Figure  5  The  explanation  for  the  poor  convergence  was  inaccurate 
approximation  of  the  implicit  operator.  The  inaccuracy  developed  from  the  combina¬ 
tion  of  the  approximate  factorization  and  the  approximate  analytical  calculation  of  the 
flux  Jacobians. 

The  standard  formulation  of  eqn(52)  allowed  the  use  of  standard  iterative  matrix 
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Figure  5:  Convergence  history  using  the  LU-SGS  method  to  invert  the  implicit  operator. 
n  is  the  number  of  physical  time  iterations,  and  m  is  the  number  of  LU-SGS  pseudo  time 
subiterations. 
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Figure  6:  Convergence  history  using  the  SGS  method  to  invert  the  implicit  operator,  n  is 
the  number  of  physical  time  iterations,  m  is  the  number  of  pseudo  time  subiterations,  and 
sgs  is  the  number  of  iterations  of  the  SGS  method. 


inversion  methods.  In  this  project  the  symmetric  Gauss-Seidel  method  was  used.  The 
convergence  history  is  shown  in  Figure  6.  n  is  the  number  of  physical  time  iterations, 
m  is  the  number  of  pseudo  time  subiterations,  and  sgs  is  the  number  of  iterations  of 
the  SGS  method.  Unlike  the  LU-SGS  method,  there  is  no  coupling  between  the  SGS 
iterations  and  the  pseudo  time  iterations.  The  values  of  the  implicit  operator  A  and  the 
inhomogeneity  vector  b  are  updated  between  pseudo  time  iterations,  but  not  between 
SGS  iterations. 

The  formulation  of  eqn(52)  into  a  standard  form  required  numerical  calculation  of 
the  flux  Jacobians.  The  limit  calculation  given  by  eqn(44)  was  simple  to  implement 
and  gave  accurate  values  of  the  flux  Jacobian.  However,  it  was  sensitive  on  the  value 
of  e  when  e  was  increased  beyond  1  x  10-12.  Using  the  complex  formulation  given  by 
eqn(46)  provided  accurate  flux  Jacobian  calculations  with  much  less  sensitivity  on  h. 
Typical  values  of  h  were  1  x  10~5. 
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Figure  7:  Mach  6  flow  impinging  on  a  hemispherical  body.  The  upper  plot  (a)  shows  the 
density  contours.  The  lower  plot  (b)  shows  the  ionization  fraction.  Notice  the  increased 
ionization  at  the  stagnation  point. 

3.3  Time  Dependent  Ionization  and  Multiple  Temperature  Ef¬ 
fects 

A  time-dependent  ionization  model  was  added  to  self-consistent ly  calculate  the  ion¬ 
ization  fraction  of  the  fluid.  The  model  is  described  by  eqn(6).  The  ionization  rate 
parameter  <  av  >ion  and  the  recombination  rate  <  av  >  recomb  were  calculated  using 
empirical  formulations  given  in  Ref.  [6]. 

The  time-dependent  ionization  model  allowed  determination  of  the  ionization  frac¬ 
tion.  For  uniform  flow  properties  the  time  dependence  is  exponential.  The  model  was 
benchmarked  against  analytical  formulations  for  its  time-dependence.  The  steady-state 
solution  was  benchmarked  against  the  Saha  equation. 

A  more  interesting  test  was  constructed  to  have  a  Mach  6  flow  impinge  on  a  hemi¬ 
spherical  body.  The  flow  was  initially  unionized,  and  ionized  upon  transition  through 
the  shock  wave.  The  results  are  shown  in  Figure  7. 

Multi-temperature  effects  have  been  added  to  the  code.  The  code  was  extended 
to  include  thermal  diffusion  for  the  constituent  temperatures  (neutrals,  ions,  and  elec¬ 
trons).  Energy  loss  mechanisms  were  added  for  the  constituent  fluid  components.  The 
constituent  temperatures  can  evolve  independently  by  diffusion  as  described  in  eqn(13). 

The  heat  conduction  equation  was  used  as  a  benchmark  to  validate  the  incorporation 
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p  =  4 


n  (b) 


Figure  8:  (a)  Strip  decomposition  and  (b)  patch  decomposition  of  a  2-D  domain. 

of  the  multi- temperature  thermal  diffusion  model. 

A  dominate  energy  loss  mechanism  for  high  electron  temperature  plasmas  is  radia¬ 
tion.  The  radiation  loss  term  due  to  Bremsstrahlung  has  been  included.  The  radiation 
was  validated  against  analytical  results. 

3.4  Parallel  Computer  Performance 

The  algorithm  was  parallelized  using  the  domain  decomposition  technique  (DDT).  The 
integral  form  of  a  general  conservation  law  was  given  by  eqn(59).  Domain  decomposition 
implementation  requires  boundary  (or  ghost)  cells  to  overlap  with  neighboring  domains 
or  blocks.  The  domain  decomposition  is  illustrated  for  two  dimensions  in  Figure  8. 

The  ghost  cells  are  used  as  boundary  conditions  to  the  real  cells  in  the  block.  A 
consequence  of  domain  decomposition  is  the  more  blocks  that  are  used  the  more  ghost 
cells  that  are  necessary.  The  ghost  cell  data  lag  the  current  computation  by  a  single 
iteration.  Therefore,  an  increase  on  ghost  cells  generated  a  slower  convergence  rate. 
Figure  9  shows  the  slightly  slower  convergence  rate.  The  convergence  history  for  4  and 
8  processors  overlay. 

The  parallel  performance  for  the  code  is  shown  in  Figures  10  and  11  for  the  code 
operating  in  explicit  and  implicit  mode.  The  grid  was  the  three-dimensional  grid  shown 
in  Figure  2  and  was  parallelized  using  domain  decomposition.  The  grid  was  scaled 
with  the  number  of  processors,  so  the  grid  size  per  processor  was  constant.  As  the 
number  of  processors  was  doubled,  the  number  of  grid  cells  was  also  doubled.  The 
ideal  speedup  was  unity.  Note  that  the  speedup  presented  is  “engineering”  speedup. 
The  value  includes  not  only  the  inefficiencies  associated  with  communication  between 
the  processors  but  also  those  associated  with  more  iterations  required  to  converge  the 
solution.  The  “engineering”  speedup  is,  therefore,  the  total  parallel  efficiency  to  obtain 
the  same  quality  of  solution  on  a  parallel  computer. 

The  high  parallel  efficiency  was  obtained  by  overlapping  communication  with  com- 
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Figure  9:  Convergence  history  using  the  SGS  method  to  invert  the  implicit  operator  on  a 
parallel  computer.  The  results  from  a  serial  computer  are  plotted  for  comparison,  n  is  the 
number  of  physical  time  iterations,  m  is  the  number  of  pseudo  time  subiterations,  and  sgs 
is  the  number  of  iterations  of  the  SGS  method. 
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Figure  10:  Parallel  speedup  for  a  three-dimensional  grid  using  domain  decomposition  on  a 
cluster  of  DEC  Alpha  workstations.  The  grid  is  scaled  with  the  number  of  processors,  so 
the  grid  size  per  processor  is  constant.  The  ideal  speedup  is  unity. 

putation.  The  boundary  information  was  exchanged  while  the  core  cells  were  computed 
and  before  the  boundary  cells  are  computed.  The  super-linear  speedup  for  the  explicit 
operating  mode  was  generated  by  slow  operation  on  a  single  processor.  This  effect  has 
been  reported  Michl  et  al  [19]  Figure  10  contains  the  performance  results  from  a  parallel 
workstation  cluster  of  16  DEC  Alpha  workstations.  Figure  11  contains  the  performance 
results  from  the  IBM  SP2  at  the  Maui  High  Performance  Computing  Center.  The 
results  on  both  computing  platforms  were  similar. 

4  Professional  Interactions 

4.1  Project  Personnel 

The  personnel  who  have  been  directly  involved  in  this  project  are  listed  below. 

_ Name _ Position _ 

Uri  Shumlak  Assistant  Professor 

D.  Scott  Eberhardt  Associate  Professor 

Thomas  R.  Jarboe  Professor 

Chris  Aberle  Graduate  Student 

John  Loverich  Graduate  Student 

Ward  Vuillemot  Graduate  Student 

Graham  Schelle  Undergraduate  Student 


20 


1  2  4  8  16 

No.  of  Processors 


Figure  11:  Parallel  speedup  for  a  three-dimensional  grid  using  domain  decomposition  on 
IBM  SP2  parallel  supercomputer.  The  grid  is  scaled  with  the  number  of  processors,  so  the 
grid  size  per  processor  is  constant.  The  ideal  speedup  is  unity. 

4.2  Collaborations 

4.2.1  Air  Force  Research  Laboratory 

Dr.  Robert  Peterkin,  Jr.  of  the  Electromagnetic  Sources  Division  of  the  Air  Force 
Research  Laboratory  at  Kirtland  AFB  on  three-dimensional  multigrid  algorithms  for 
MACH3,  development  of  a  parallel  PIC  (particle  in  cell)  code  for  microwave  simula¬ 
tions,  and  stabilization  of  the  Rayleigh-Taylor  instability  in  solid  liner  implosions  by 
introducing  a  sheared  axial  flow  by  designing  conical  liners.  Knowledge  developed  in 
the  area  of  relaxation  schemes  was  implemented  into  the  ICEPIC  code  to  make  a  3-D 
Poisson  solver.  The  solver  was  needed  to  determine  electric  field  concentration  on  a  high 
power  microwave  source.  The  collaboration  occured  in  person  during  August.  Several 
phone  and  email  discussions  took  place  throughout  the  year. 

4.2.2  Sandia  National  Laboratories 

Dr.  Norm  Roderick  of  the  Pulsed  Power  Sciences  Center  at  Sandia  National  Laboratories 
on  the  uses  of  sheared  axial  flows  to  stabilize  z-pinch  implosions.  This  is  an  ongoing 
collaboration  that  resulted  in  the  publication  listed  in  the  following  section. 

4.2.3  University  of  Washington 

Prof.  Scott  Eberhardt  of  the  Aeronautics  and  Astronautics  Department  and  Prof.  Randy 
LeVeque  of  the  Applied  Math  Department  on  approximate  Riemann  solvers  and  their 
applications  to  multidimensional  problems.  We  have  regular  discussions  on  a  weekly 
basis. 
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Prof.  Tom  Jarboe  of  the  Aeronautics  and  Astronautics  Department  on  the  higher 
mode  stability  of  spheromaks  and  on  the  effect  of  realistic  three-dimensional  geometries 
on  spheromak  stability.  This  is  an  ongoing  collaboration  that  resulted  in  the  publica¬ 
tions  listed  in  the  following  section. 

4.3  Publications 

A  journal  article  describing  our  algorithm  has  been  published  in  the  Journal  of  Compu¬ 
tational  Physics .  The  title  is  “An  Implicit  Scheme  for  Nonideal  Magnetohydrodynam¬ 
ics”  by  O.  S.  Jones,  U.  Shumlak,  and  D.  S.  Eberhardt.[13]  The  citation  is  Journal  of 
Computational  Physics  130,  231  (1997).  Another  journal  article  which  describes  the 
use  of  sheared  flows  to  stabilize  the  Rayleigh-Taylor  instability  has  been  published  in 
Physics  of  Plasmas.  This  is  work  that  was  performed  with  collaboration  at  the  Air 
Force  Research  Laboratory.  The  title  is  “Mitigation  of  the  Rayleigh-Taylor  Instability 
by  Sheared  Axial  Flows”  by  U.  Shumlak  and  N.  F.  Roderick,  [20]  The  citation  is  Physics 
of  Plasmas  5,  2384  (1998). 

Two  papers  describing  the  higher  mode  stability  in  spheromak  plasmas  and  the 
effect  of  realistic  three-dimensional  geometries  on  spheromak  stability  have  also  been 
published. [21,  22]  The  citations  are  Physics  of  Plasmas  6,  4382  (1999)  and  Physics  of 
Plasmas  7,  2959  (2000). 

Two  papers  describing  this  project  will  be  presented  at  the  upcoming  AIAA  Com¬ 
putational  Fluid  Dynamics  conference.  One  paper  titled  “An  Approximate  Riemann 
Solver  for  MHD  Computations  on  Parallel  Architectures”  will  present  an  overview  of 
the  project,  and  the  other  paper  will  present  recent  work  on  the  analytical  flux  Jacobian 
calculation. 


5  Conclusions 

The  successful  development  of  the  three-dimensional  advanced  implicit  algorithm  and 
the  implementation  of  time-dependent  ionization  and  multiple  temperature  effects  show 
that  this  project  is  reaching  its  objectives.  The  research  related  to  this  project  has  been 
published  in  refereed  journals  and  presented  at  international  conferences.  Valuable  col¬ 
laborations  have  been  formed  with  the  Air  Force  Research  Laboratory,  Sandia  National 
Laboratory,  and  other  universities. 

The  continuing  development  of  this  project  will  include  investigating  the  TVD  prop¬ 
erties  and  the  analytical  flux  Jacobian  methods.  An  important  result  of  this  work  is 
the  determining  of  the  difficulty  in  stabilizing  the  Hall  effect.  A  two  fluid  plasma  model 
may  be  essential  to  properly  treat  Hall  effect  plasma  physics.  [23] 
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